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A unique numerical method has been developed for solving one-dimensional ablation 
heat transfer problems. This paper provides a comprehensive description of the method, 
along with detailed derivations of the governing equations. This methodology supports 
solutions for traditional ablation modeling including such effects as heat transfer, material 
decomposition, pyrolysis gas permeation and heat exchange, and thermochemical surface 
erosion. The numerical scheme utilizes a control-volume approach with a variable grid to 
account for surface movement. This method directly supports implementation of non- 
traditional models such as material swelling and mechanical erosion, extending capabilities 
for modeling complex ablation phenomena. Verifications of the numerical implementation 
are provided using analytical solutions, code comparisons, and the method of manufactured 
solutions. These verifications are used to demonstrate solution accuracy and proper error 
convergence rates. A simple demonstration of a mechanical erosion (spallation) model is also 
provided to illustrate the unique capabilities of the method. 

Nomenclature 

a = nodal coefficient 

A = area, m 2 ; Arrhenius pre-exponential, sec" 1 ; slope of linear in-depth flux versus surface temperature 
b = nodal coefficient 

B = offset of linear in-depth flux versus surface temperature 

B' = nondimensional mass flux of ablation products away from the surface 

B[ : = nondimensional mass flux of pyrolysis products into the surface 

S c = nondimensional mass flux of surface material erosion 
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c 


nodal coefficient 


c p = specific heat at constant pressure, J/kg- K 
C = generic capacitance coefficient 

C H = Stanton number for heat transfer 

Cm = Stanton number for mass transfer 

d = nodal coefficient 

E = activation energy, J/kmol 

F , = curvature factor, m 

h = heat transfer coefficient W/nr-K; sensible specific enthalpy, J/kg 

h = effective specific enthalpy, J/kg 

H r = recovery enthalpy, J/kg 

I e = Arrhenius integral 

j k = diffusion mass flux of the k th element summed over all species, kg/m 2 -sec 
K = mass fraction; generic diffusion coefficient 

K k = mass fraction of the k' h element summed over all species 

k = thermal conductivity, W/m-K 

Le = Lewis number 

m = reaction order 

m = mass flow rate, kg/sec 

m" = mass flux, kg/m 2 -sec 

m = volumetric mass generation, kg/m’-scc 
MW = molecular weight, g/mol 

N = the node number corresponding to the back surface 

N SU rf = the node number corresponding to the front (ablating) surface 

P = pore pressure. Pa 

q" = heat flux, W/m 2 

(f = volumetric energy generation rate, W/m 3 



Q = energy generation rate, J/sec 

Q p = heat-of-pyrolysis, J/kg 

r = radius, m 

R = ideal gas constant, J/g-K; nondimensional nodal coefficient modifier for radius correction 

R„ = universal ideal gas constant, J/mol-K 

s = surface location, m 

s = erosion rate, m/sec 

S = generic source term 

Sc = constant portion of a generic source term 

Sp = linear coefficient in generic source term 

t = time, sec 

T = temperature, K 

T c hem = combined parameter used for surface heat flux calculations, J/kg 
u = velocity of boundary layer gas, m/sec 

v D = Darcy velocity, m/sec 

V = velocity of pyrolysis gases, m/sec 

x = depth from the original location of the front surface, m; fraction of virgin mass loss associated with the 

; th pyrolyzing component 

Z* = mass transfer driving potential 

Greek 

a = overall extent-of-reaction 

or, = extent-of-reaction for ; th pyrolyzing component 

8x = nodal spacing, m 

Ax = cell thickness, m 

At = time step, sec 

r = permeability, m 2 ; thermal diffusion coefficient, kg/m-s 

£ = modified non-dimensional ablation rate used in blowing corrections 


K 


total number of chemical elements 



A = parameter used in blowing correction model 
/j = dynamic viscosity, N-sec/m 2 

p = bulk density, kg/m 3 

p = actual density, kg/m' 

pc = effective capacitance, J/m 3 -K 


(f> = volumetric porosity 

ip = generic solution variable 


Subscripts 

adv = 
b 

c = 
cond = 
conv = 
chem = 

D 

e = 

f = 

fe-g- = 
g 

gen = 
i = 
inc = 
k = 
M = 

N 

o = 
old = 


advection 

back surface 

fully-charred 

conduction 

convection 

thermochemical 

related to Darcy’s model 

boundary layer edge; on the “east” side of the cell or node (away from the front surface) 
formation 

with chemical composition at the boundary layer edge (frozen edge gas) 

pyrolysis gas 

generation 

the ; th pyrolyzing component; the ; th specie 
incident 

the k' b chemical element 
boundary layer mass transfer 
last element number 
conditions in the absence of blowing 
at the previous time step 



p 


constant pressure 


P = the local node or cell 

r = recovery 

ref = reference value 

s = solid material; the solid surface 

sto = storage 

surf = surface 

tch = thermochemical 

v = fully-virgin material 

w = at the surface (wall); on the “west” side of the cell or node (toward the front surface) 

4w = the four-wall heating model 

i// = generic solution variable 



Introduction 


Ablative insulators are commonly used in aerospace applications to protect structural components from extreme 
aerothermochemical environments. The physical phenomena associated with ablation heat transfer depend on the 
application, but most involve in-depth material pyrolysis (charring) and thermochemical surface ablation. Related 
physics are summarized in the pioneering work done by Aerotherm in the 1960s [1] - [6]. Numerical modeling 
requires the solution of an energy equation including the effects of pyrolysis on a domain that changes as the surface 
ablates. In addition, surface movement may occur due to material spallation [7], [8] and intumescence [9]. These 
effects further complicate the numerical treatment of the changing domain. A numerical modeling program (ITRAC) 
has been developed that includes capabilities for modeling these and other complex ablation phenomena [10] - [12]. 
At the foundation of the program is a variable grid control-volume method for numerical solution of the governing 
equations. This paper provides details on the mathematical models governing ablation phenomena along with a 
description of the numerical solution methods used in the program. In addition, verification cases are presented that 
support accuracy of the method and proper program implementation. 

Various methods exist for the solution of moving boundary heat transfer problems, often referred to as “Stefan 
problems” after the early work of Stefan [13]. Crank [14] classifies related numerical methods into two types; 
“front-tracking” and “front-fixing.” With front-tracking methods, the ablating surface (front) is tracked as it moves 
into the material, and the spatial discretization is updated in some manner to account for the changing domain. This 
may be in the form of a completely new grid or a modification of the mesh near the ablating front. Crank describes 
difficulties using front-tracking methods associated with requirements for equal grid spacing in finite-difference 
schemes. He discusses modified finite -difference methods and special temporal discretization schemes that 
overcome these difficulties. With front-fixing methods, the location of the moving boundary is fixed in the solution 
domain by using a transformed spatial coordinate such as that of the commonly used method of Landau [15]. With 
this method, the ablating surface is fixed in the transformed coordinate and the spatial grid is unchanged at each time 
step. The influence of the moving boundary is accounted for in an additional advective term that makes its way into 
the transformed energy equation. Use of these methods can result in awkward interfacing between ablating materials 
and nonablating substrate materials. In addition, modeling of instantaneous spallation is not possible since the 
advective term in the transformed energy equation becomes infinite. 



The most commonly used program for ablation modeling in aerospace applications is CMA [16], which was 
developed as part of the work by Aerotherm. CMA relies on surface thermochemistry tables generated by the 
companion ACE program [17]. Discretization in CMA is based on a transformed energy equation for the ablating 
material with special treatment of a transition node between ablating and nonablating materials. The transition node 
changes in size to accommodate surface recession and is dropped when the size reaches a specified minimum. More 
recently, Amar et al. [18] used a front-fixing method with a Landau transformation to successfully solve pyrolysis 
and thermochemical ablation problems. Their method reduces complexity by using a contracting grid that does not 
require nodal dropping at the interface. The work of Amar et al. also provides some consideration of program 
verification, which has been sparse within the available ablation modeling literature. 

This paper gives details of the fundamental equations governing ablation heat transfer phenomena. The 
derivations include descriptions of related modeling assumptions and property definitions. In addition, details are 
provided on the general numerical solution method used to solve these coupled equations. As opposed to the more 
commonly used front-fixing methods, the approach here is based on a front-tracking scheme with a variable grid. 
The method simplifies treatment at material interfaces and does not require special discretization at an interface with 
a nonablating material. Nonuniformity in nodal spacing is naturally handled using the control-volume method of 
Patankar [19], which does not require uniform nodal spacing. The method allows simple implementation of complex 
erosion models, such as spallation, since nodes (and control-volumes) are removed from the front surface in a 
manner that mimics the actual physical process. Material swelling is also easily implemented. The ITRAC program 
supports various models for thermochemical and mechanical ablation, in-depth material pyrolysis, and material 
swelling. One -dimensional planar, cylindrical, and spherical coordinates are supported. Complete description of all 
the solution features in the program is beyond the scope of this paper, which focuses on the mathematical models 
and numerical solution schemes supporting general ablation heat transfer problems. The numerical approach is 
described in detail, and solution accuracy is verified using analytical solutions, manufactured solutions, and 
comparisons with accepted codes. While the resulting verification does not provide exhaustive assessment of the 
many features of the program, it does provide verification of the core solutions. 



Phenomena 


The primary phenomena associated with ablative insulators are illustrated in Fig. 1. The insulator can be heated 
by radiation and convection heat transfer at the front surface. As regions within the insulator increase in temperature 
the material decomposes (pyrolyzes), and a pyrolysis front progresses into the insulator leaving behind a layer of 
charred material. It is customary to define pyrolysis and char depths corresponding to the onset and completion of 
pyrolysis; the region between the two depths is then referred to as the pyrolysis zone. Pyrolysis gases are generated 
as the material chars. These gases flow through and exchange energy with the porous char structure. Meanwhile, 
erosion of the surface material can occur due to chemical and mechanical interaction at the adjacent boundary. 
These primary phenomena are not inclusive of all of the complex physics that can generally influence ablation heat 
transfer. Some key assumptions and simplifications in the present modeling are: 1) internal gas generation is only 
due to material pyrolysis, 2) pyrolysis gas transport is driven only by permeation, 3) pyrolysis gas enthalpy is a 
unique function of temperature, 4) local thermal equilibrium exists between the pyrolysis gases and the porous char 
material, 5) thermal and chemical equilibrium prevail at the material surface, and 6) intermediately charred material 
properties are uniquely defined by the degree of char. As a result, the described model does not include such effects 
as gas generation from internal volatiles, internal multi-species transport and related chemical reactions, kinetically- 
controlled reactions within the pyrolysis gases, the effects of in-depth condensation of pyrolysis products (coking), 
or the influence of heating rate on material properties. The significance of these assumptions and simplifications are 
jj jcific to different applications and materials. However, the primary phenomena identified here are common to a 
Side range of ablation heat transfer applications, and form a foundation upon which additional fidelity can be added. 

Mathematical models are described below for these transient thermal, pyrolysis, and ablation phenomena. The 
models are for one-dimensional planar, cylindrical, or spherical geometries. Related solution variables are the 
temperature field T(x,t), the extent-of-reaction (degree-of-char) field a(x,t). the pore pressure field P(x,t), and the 
location of the ablating surface s(t). In the present formulation, the eroding surface is labeled as the “front” surface 
and the “back” surface is fixed as shown in Fig. 2. Details of the mathematical models along with their numerical 
solution are presented in the sections below. 


Various uses can be found in the literature for the terms “erosion” and “ablation.” Here “erosion” is used to refer to the loss of surface material 
due to chemical or mechanical effects, and “ablation” refers to total mass loss including the effects of both surface erosion and in-depth pyrolysis. 


Governing Equations 


Presentation of the governing equations begins with a description of models for material pyrolysis (charring) that 
occurs as the material is heated. This process is modeled using a bulk solid density p s that decreases from an initial 
“virgin” value p v to a final “fully charred” value p c as pyrolysis progresses. The extent of this conversion process is 
quantified using an extent-of-reaction a that progresses from an initial virgin value of zero to a fully-charred value 
of 1. 

The second section describes models for pyrolysis gas permeation through the porous char structure generated in 
the charring process. These models capture the important physics of pyrolysis gas flow that affects in-depth thermal 
transport as well as thermal and chemical interactions at the material surface. The related governing equations are 
based on mass conservation applied to the gas phase. The primary solution variable is the pressure field P within the 
porous material. 

The third section describes energy transport within the material. Mechanisms include thermal conduction and 
storage, generation effects of material conversion (charring), and heat exchange between the porous solid and the 
permeating gases. The related solution variable is the temperature field T within the material. 

The final two sections describe surface phenomena related to heat transfer and thermochemical erosion. Related 
models in these sections provide important thermal boundary conditions as well as models for the final solution 
variable s representing the location of the ablating surface. Details are given below. 


Pyrolysis Kinetics 

The extent of material pyrolysis is quantified using an overall extent-of-reaction a based on the bulk density of 
the decomposing solid p s related to densities in the fully-virgin and fully-charred conditions p v and p c . The 
relationship is 


a = 


Py-Ps 


(1) 


Pv - P, 

It is convenient in the derivation of the governing equations to express the rate of change of solid density in 
terms of the rate of change of a. From Eq. (1) the relationship is 


= -(Py-P C ) 


da 

dt 


OPv 

dt 


(2) 



The pyrolysis process is modeled as a combination of multiple reactions, each with its own extent-of-reaction a-,. 
The fraction of the total virgin mass loss as a result of the ; th reaction is denoted x,. Assuming that the material 
volume is constant during the decomposition process, the bulk density of the decomposing solid is related to the 
component reactions by 

Ps= A.- PvYs x i a i ( 3 ) 


A fully-charred condition corresponds to a value of unity for each a, giving the following for the fully-charred 
density 

Pc=Pv~PvT. X i ( 4 ) 


Combination of Eqs. (1), (3), and (4) gives the following expression for the overall extent-of-reaction 





( 5 ) 


The rate for each component is modeled using the following Arrhenius expression 

8c ^=A,e- E ‘ ,RT (l-a,p 
dt ' 

According to Eq. (5), the overall rate of change is related to the individual rates by 


Yr 

da Y ‘ dt 


dt 




(6) 


( 7 ) 


In-Depth Mass Balance 

A control-volume for gaseous mass balance within the porous pyrolyzing material is shown in Fig. 3. The figure 
shows a gaseous mass generation term m’” en associated with pyrolysis, advection of mass in the form of the mass 

flux m " g , and storage of gas having density p g in the solid pores. Considerations of the storage, advection, and 

generation terms are described below. The resulting terms are combined into a governing mass balance equation. 
Storage 


For a cross-sectional area A the storage rate is 



where (p is the local volumetric porosity. 

Advection 

The local mass flux m " is defined as positive when flowing toward the surface against the positive x-direction. 
In terms of the mass flux, the net advection into the control-volume is then 


m adv,net 


dx 


(9) 


Generation 


The volumetric gas generation rate m gen is equal to the negative of the rate of change of solid density dp s / d t . 
This, along with Eq. (2), gives 


dp 


da 


m, e n = m”Adx = — -p-Adx = (p v - p c )—Adx 


dt 


dt 


( 10 ) 


Mass Conservation Equation 

Conservation of mass for the gaseous phase requires a balance of storage, net advection, and generation such that 


™ sto ™adv,net ™gen 


(11) 


Combination of Eqs. (8) - (1 1) gives the mass balance as 


d(P g P) 

dt 


1 d(rit’A) 
A dx 


+ (Pv 



Ideal gas behavior is assumed so that 


P e 


P 

~RT 


( 12 ) 


(13) 


In terms of the Darcy velocity v D , also defined as positive against the x-direction, the mass flux can be written as 

m" g = P g v D (14) 

According to the Darcy model, the Darcy velocity can be written in terms of the pressure gradient as follows 


J ^_dP_ 

p g dx 


(15) 


where T is the material permeability, //,, is the pyrolysis gas viscosity, and the positive sign on the right-hand-side 



of the equation is due to the directional definition of v D as positive against the x-direction. 


Porosity is modeled as a linear function of the degree-of-char as follows 

<t = fa(\-a) + <!> c a 

Incorporating Eqs. (13) - (16) into Eq. (12) gives the following for the gaseous mass balance 


<j> dP _ 1 8 
RT dt A dx 


p r 8P 

RT ju g dx j 


P<j> dT 
RT 2 dt 


+ Pv~Pc + 


p<t> v P<t>c 1 da 


RT RT ) dt 


For planar, cylindrical, and spherical coordinates this becomes 


<f> dP _ d 
RT dt dx 


<t> dP _ 1 d 
RT dt r dx 


dP 1 d 


P T 8P 
RT /j dx 




p r dp 

RT ju dx 




Ptj> dT 
RT 1 dt 
planar 

P<j> dT 


+ Pv~ Pc + 


p<t> v P<t>, \ da 


RT RT J dt 


RT 2 dt 
cylindrical 


+ Pv -Pc + 


p<t> r P<t> c ^ da 


RT RT j dt 


RT dt 


dx 




P T dP 

RT p g dx j 


P<j> dT 
RT 1 dt 
spherical 


+ Pv - Pc + 


P<t>v P<t>c 1 da 


RT RT ) dt 


(16) 


(17) 


(18) 


(19) 


( 20 ) 


The radius r used in Eqs. (19) and (20) is related to the local x-position according to the concavity and the radius of 
the front surface r s . For a concave geometry, the radius increases with x, for a convex geometry it decreases. The 
corresponding relations are 

(21) 


r = r s + x 

concave geometry 

r = r s -x 
convex geometry 


(22) 


In-Depth Energy Balance 

A control-volume for energy balance within the porous pyrolyzing material is shown in Fig. 4. Energy is stored 
in both the solid and gas phases within the control-volume, while energy enters and leaves through gaseous 
advection and conduction. In addition, the pyrolysis process contributes volumetric generation (usually 
endothermic) as illustrated. Each contribution to the balance of energy within the differential element is considered 


below. 



Storage 


For a cross-sectional area A the storage rate is 


,, d. , „ , „ ,dh d(p</>) dh, , dp , 

a,. = + A". >-"* = ".'ir + *« -ft- + ~dt* h, ~dt 


+ p — - + /7 — - Adx 
s dt s dt 


Conduction 


Net conduction into the control-volume is 


a d iA dT \n 

Qcond= T" kA—\dX 
C/X I ux ) 


Net advection into the control-volume is given by 


d(m" Ah ) T „ dli d(m"A ) 

Qadv = dx = ^l-A — + h„ dx 

ad ' dx 8 dx g dx 


Generation 

Generation is a function of the pyrolysis rate and the heat-of-pyrolysis Q p . With Q p a positive value for 
endothermic pyrolysis, the generation rate is 


Q gen = <t'"en Adx = Q P ~^~ Adx 


Energy’ Conservation Equation 


Conservation of energy requires the following balance of storage, conduction, advection, and generation 


Qsto Qcond Qadv Qgi 


Combination of Eqs. (23) - (27) gives 


dh g d{p </>) dh s . dp. Id (..dT 
dt g dt dt dt A 5x1 dx 


- + p s ^ + h s ^ = -— kA— 
dt dt A dx l dx 


. „ oh h d(m"A) dp, 

+m " — + — — + Q n -oo- 

8 dx A dx p dt 


After combining with Eqs. (2) and (12) this can be written as 


« dh dh 1 d ( dT ) dh doc 

p 6 — — + p — — = kA — + rh — — — {Q —h +h )(p —p ) — 

PgV dt p ‘ dt Adx { dx) % dx s gAPv Pc> dt 


The solid enthalpy /?, is assumed to be a linear function of the overall extent-of-reaction a based on fully-virgin 


and fully-charred values h v and h c as follows 



h s =h v (\-a) + h c a 


(30) 


This gives the following for the rate of change of solid enthalpy 


dh s _ dh s dT dh s da 
dt dT dt da dt 


dh . dh c 

— L (1 - a) h -a 

dT dT 


dT , ,6c( 

— + (-h+h.) 

dt ' c dt 


Combining Eq. (31) with Eq. (29) and incorporating specific heats ( c p =dh/dT) gives 
[Pgfcp.g + Ps( ] ~ a ) C p,v + Ps^p.c \^~ 


A dT\ 

.. dh 

Q„ - 

f h + p s(K-K)} 

\ + h 

■ kA — \ + 

\ dx) 

8 dx 


V Pv Pc J 

S 


(A - Pc) 


da 

dt 


where thermal equilibrium has been assumed between the gas and solid phases. Effective capacitance and 
enthalpy terms are defined as follows 

pc = Pg^p.g + Ps (1 - a )c p , v + P s ac p ,c 

and 

, 7 .* 

Pv - Pc 

Incorporating these into Eq. (32) gives the following general equation representing energy conservation 


pc 


dt Adx{""dxJ "' 8 dx ( ' Qp k +hg)(Pv Pc) - 


. da 
dt 


For planar, cylindrical, and spherical coordinates this becomes 

— dT d(,dT\ . „ dh 


pc— = —I k— | + m " — — - (Q -h + h )(p v - 
dt dx v dx J 8 dx ‘ g c dt 

planar 

— dT 1 d ( dT\ . „ dh - da 

^ = 7& (’* *> "* a7 “ ,e ' “ " + *' )(ft - »77 

cylindrical 

— dT 1 d( 2, dT'] ,„dh g a 

spherical 

with radius definition given in Eqs. (21) and (22). 


(31) 

(32) 
solid 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 


Surface Energy Balance 



Kendall et al. [20] provide a model for coupling chemically reacting boundary flow to an ablating material. This 
approach is predominantly used within the ablation modeling community. With this model, coupling is based on 
transfer coefficients to capture the effects of species diffusion through the concentration boundary layers and heat 
transfer through the temperature boundary layer. Surface conditions, including the surface heat flux and ablation 
rate, are determined based on the assumption of chemical equilibrium at the surface temperature and pressure. The 
surface temperature, an unknown in the modeling, must be determined in a manner that reconciles surface energy 
conditions with an in-depth energy balance. Three modeling options are provided based on different conditions 
within the reacting boundary gasses. The options correspond to conditions of 1) unity Lewis number Le, 2) non- 
unity Le, but equal diffusion coefficients for the various species within the boundary flow, and 3) nonunity Le and 
unequal diffusion coefficients. The models are based on the energy terms depicted in Fig. 5. Brief descriptions are 
given in the sections below which categorize the terms according to thermochemical convection, radiation, mass 
transfer, and in-depth conduction. A final subsection discusses the overall surface energy balance equations. 
Thermochemical Convection 

Heat transfer from a chemically reacting flow includes the effects of energy transfer through a temperature 
gradient and enthalpy transfer from species diffusion through concentration gradients. Within the ablation modeling 
community, these separate forms of energy flux are often referred to as "convective" and "chemical" heat fluxes. 
Here, the combined effects are referred to as "thermochemical convection," and the term q" ch is used to denote this 
type of energy transfer. Different transfer coefficient models are used according to the assumed conditions in the 


boundary gases. For the case of unity Le, thermochemical convection is modeled by 


(39) 


%:h = Pe U e C H ( H r ~K) 


unity Le 

For nonunity Le with equal diffusion coefficients the model becomes 


‘ill, = Pe U eC H (H r ~ h w ) feg + P e U e C M Yj K ie ~ Z K m ] A /* 

i i J 

equal diffusion coefficients 


For the general case of nonunity Le and nonequal diffusion coefficients the model is 


(40) 


‘Itch ~ Pe U e^ r ^w) fe.g + P e U M ^ ^ ie ^ ^ h 

Vi i J 

equal diffusion coefficients 


h 


(41) 


Radiation 


Radiation from adjacent boundary flow is common in ablation heat transfer applications. Here, incident radiation 
absorbed at the surface or "wall" is denoted by oc w q" md inc , with a v representing the wall absorptivity. The surface 

also reradiates according to the wall emissivity £ w and temperature T w . Hence the term s w oT 4 in Fig. 5 to account 
for reradiation. 

Mass Transfer 

Energy is transferred with mass transfer to and from the eroding surface. Mass flow into the surface includes the 
mass flux of pyrolysis gases ih " and the mass flux of charred surface material m" consumed in the surface 


reactions. The gross mass flux of reaction products away from the wall is given by ( pV) w . Related enthalpies are 


lose of the pyrolysis gases h g , the charred surface material h c , and the wall reaction products h w . Energy fluxes 


into the surface are then m" g h g and m"h c , while energy out of the surface is (pV) w h H , 


In-Depth Conduction 


The energy conducted into the solid material is given by q” cond 


Energy’ Balance 

Ihe surface energy terms are combined into an overall balance given by 

q",C + «w q"ra d ,in C ~ S « (jT ' ~ ( P V )w K + >K k c + ™" g h g “ = 0 ( 42 ) 

Incorporating the models of Eqs. (39) - (41), and rewriting the balance with q” ond on the left-hand-side of the 


equation gives: 


(43) 


q" d = puCJH -a + B')h +B'h +B'h ] + a q" -s aT 4 

™cond ree Hi- r V ' w c c g g-l w “rad,inc w w 

unity Le 


q'cond =Pe U e C ,A H , 

~K)f, 

g 




+Pc U e C M 

I?'* 

“Z K^ll" + B'A + B' g h s - B'h w 

+ «w <l"radM 


(44) 


equal diffusion coefficients 


q'cond ~ P e U e G h(H , 

hf) fe. g 


+P e u e C M 

(Z Z -> -ZCV +B 'K + B 'ghg - B'h w 

„ r4 (45) 

+ «w q rad ,mc - £ y, (jT ^ 

unequal diffusion coefficients 








The models are written in terms of nondimensional mass fluxes of pyrolysis gas, surface (char) consumption, 
and total gaseous flow away from the wall ( B [, , B , and B ' ). The nondimensional “B-prime” definitions are given 
below 


m 

ZV S 

s — n 

Pe U e C M 


B' =- 


m„ 




B'= (pV ' 1 " =B' +B[. 

Pe U e^M 


(46) 

(47) 

(48) 


The parameter B is commonly used to represent a nondimensional blowing rate at the wall [21]. The non- 
dimensionalization is usually based on mass transfer coefficients in the absence of blowing. The prime notation in 
the definitions of Eqs. (46) - (48) is used to denote nondimensionalization based on mass transfer coefficients that 
have been adjusted for blowing effects. The adjustment is based on correlations for transpiring boundary layers and 
is given by the following relation [16] 


C H _ j 
C h„ e c -l 


(49) 


with 


2 Am" 

Pe U e^H 0 


(50) 


and 


• ft • ff , • ft / c 1 \ 

m =m c +m g (51) 

In these equations, C H is the Stanton number, C Ho is the Stanton number prior to the blowing adjustment, that is, the 
value undisturbed by mass injection, m " is the total mass flux of gaseous products injected into the boundary layer, 
and A is a parameter of the semi-empirical model used to account for different flow conditions. A value of 0.4 is 
commonly used for turbulent flows; with a value of A = 0.5, the correlation reduces to that provided by Kays and 


Crawford [22]. 



The various parameters required in Eqs. (43) - (45) are provided in thermochemical (B-prime) tables that list the 
parameters as functions of El„ and S c . These tables are generated, for example, using the ACE surface 
thermochemistry program [17]. 

Surface Erosion 

The selected surface energy balance model of Eqs. (43) - (45) must be reconciled with the in-depth energy 
solution. Once this is completed (at each time step), the thermochemical erosion rate is calculated from the resulting 
B' c as follows 

■ = KPeU e C M (52 ) 

Pc 

Numerical Solutions 


Domain Discretization 

The domain is discretized into control-volumes following the method of Patankar [19]. An illustration is 
provided in Fig. 6. Each control-volume is labeled as CV t and corresponding nodes are labeled as N t . Interfaces 
between control-volumes are denoted by /,. The initial grid is defined by specifying the thickness Ax t of each 
control-volume. Corresponding nodes are centered within the control-volumes with the exception of the front and 
back boundary nodes, which are placed at the edge as shown. Nodal spacings Sx t are calculated accordingly. 
Numbering of the control-volumes and nodes begins at t = 1 for the front surface and increases in-depth to a final 
value of B for the back boundary. 

Variable Grid Methodology 

The grid is modified for surface ablation as illustrated in Fig. 7. The location s of the ablating surface for each 
time step is determined based on the calculated ablation rate s from the previous time step. As time progresses, s 
moves deeper into the domain and the grid is partitioned into active and inactive parts as shown. Solutions are then 
performed only on the active portion of the grid, which continuously decreases in size as control-volumes are 
deactivated. At each time step, an assessment is made as to which control-volume contains the new surface. Once 
this is determined, control-volumes at the surface are adjusted and the new surface node is labeled N s - This defines 
the active portion of the grid to be associated with nodes N s through N B . Adjustments of the surface control-volumes 


are made according to where s falls within the discretized domain. Four cases are considered as shown in Fig. 7, 
which shows the old grid on the left and the modified grid on the right. Each case is discussed below. 

Case I 

Flere the surface falls within CV, and above TV,- (Fig. 7a). CV t is split as shown on the right of the figure. TV,- is not 
moved. N s is identified as TV,-_i which is repositioned at the new surface location ( x = s). The upper boundary of CV t . i 
is defined by the new surface location. The interface / M is repositioned midway between TV,-. and TV,-. Ax h dx f _ i and 
Ax /. i are modified as shown. 

Case II 

Flere the surface falls directly on TV,- (Fig. 7b). In this case TVs is identified as TV,, and the lower half of the old CV 
becomes the new CV,. The location of TV,- is unchanged so that it is now positioned at the upper boundary of CV,. The 
location of I, is unchanged. Ax, is adjusted accordingly. 

Case III 

Here the surface falls below TV, (Fig. 7c). N s is identified as TV,-, and the lower portion of the old CV, becomes the 
new CV h Nj is moved to the upper boundary of CV, (the new surface location), and the location of I, is unchanged. 
Ax, and dx, are adjusted accordingly. 

Case IV 

Here the surface falls directly on /,- (Fig. 7d). N s is identified as TV,-, and CV, is split into the new CV, and CV,-/ as 
shown. TV,- is moved to the upper boundary of the new CV,- (the surface location), and /, is repositioned midway 
between TV,- and TV, + 7 . Ax t and Sx t are adjusted accordingly. 


Extent-of-Reaction Solutions 

Extent-of-reaction solutions depend only on the local temperature history (Eq. ( 6 )). As a result, local solutions 
can be obtained at each node by direct integration of Eq. ( 6 ) with respect to time. This is done at each time step after 
separating variables. The equation for the z th extent-of-reaction a, is then 


P (\-a i y m ‘da i =A i \‘e~ El,RT dT 

J a, „ J t„ 


(53) 


where a i 0 is the “old” value (at the beginning of the time step) and a, and r are dummy integration variables. The 
term on the right-hand-side of Eq. (53) is the Arrhenius integral I E i , that is 




Numerical difficulties can be associated with Eq. (56). In particular, if the reaction order m, is less than 1.0, 
inaccuracies and stability issues can arise if time steps are too large. Consider the special case of m, = 0.5, for which 
the exponent on the bracketed term of Eq. (56) has a value of 2.0. This case is illustrated in Fig. 8, which 
corresponds to a i0 = 0.0, a particular value for I Ei , and various time steps. For these conditions, the results become 
unrealistic for time steps beyond 0. 1 sec since a, values begin to reduce, even becoming negative, with increased 
time steps beyond that value. The term in brackets from Eq. (56) is also shown in the figure. It can be seen that the 
region beyond the peak a, value corresponds to negative values of this quantity. For other reaction orders (< 1.0), 
the bracketed term becomes complex, resulting in stability problems. In order to correct for these conditions, the 
bracketed term is monitored and if it becomes negative, the calculation is bypassed and the maximum a , value of 1 .0 
is applied. 

Numerical difficulties can also be encountered for m, = 1.0. If I Ei becomes too large, the exponential term in Eq. 
(57) can exceed machine precision, resulting in an “infinite” result. In order to avoid this, an upper limit of 50.0 is 


imposed on I E i . 



Pore Pressure and Temperature Solutions 

The general governing equations for pressure and temperature are the mass balance of Eqs. (18) - (20) and the 
energy balance of Eqs. (36) - (38). These equations can be summarized in terms of a storage term, a diffusion term, 
and a source term as follows 


Cr-^ = 4-(r-Kf)HSc + S„ry 


(59) 


dt 8x V 

where n has the value 0, 1, or 2 for planar, cylindrical, and spherical conditions respectively, and the source term is 
written as a linear function of the solution variable. Table 1 gives corresponding definitions for (//, C, K, Sc, and S P 
for the mass equation and the energy equation. Discretization equations are derived for each node (cell) of the 


domain. Special considerations are made for internal nodes and the boundary nodes. In the discretization of the 


governing equations, the solution variable is generally treated as temporally uniform and spatially linear. This results 
in a scheme that is first-order in time and second-order in space. Details are given in the sections below. 

Internal Nodes 

Discretization equations are generated for each internal node by integrating the governing equation, Eq. (59), in 
space (over the corresponding cell) and time (over the time step). Special considerations are made based on 
concavity, which is illustrated in Fig. 9 for a cell surrounding the P th node; cell boundaries are denoted by the dotted 
lines. The integration is given by 

re r t+At 8\I/ re r t+At 8 1 8w\ r' f <+Af r « f >+St 

ff Cr" ^-dtdx = f f —\r n K^-\dtdx+ ff S r r"dtdx+ ff S p wr"dtdx (60) 

JvrJr Qt JwJf <3 r l Q r J JwJ t c J , vh pr v 7 


Each of the four integrals is considered separately below. 

The first integral in Eq. (60) represents capacitance. Here, temporal integration is based on the assumption that C 
is fixed over the time step. In addition, the solution value is assumed to be uniform over the cell. With these 
assumptions, the integral is 

j, ' r 4f Cr " D g t dtdx = CAx p F r (¥ P -Y° P ) (6 1 ) 

Where F r is a correction factor that accounts for curvature and has the definitions listed in Table 2. 

The second integral in Eq. (60) represents diffusion. Here, temporal integration is based on the assumption that K 
is fixed over the time step. The solution variable is treated implicitly, with new values prevailing over the time step. 


Finite difference approximations are used for spatial derivatives of the solution variable. The resulting integral is 


' dr 
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(62) 


8x e 8x„ 

where K e and K w are values of K at the interfaces of the cell. Following Pantakar, these values are calculated as the 
harmonic mean of values evaluated at the nodes on either side of the interface. 

The third integral of Eq. (60) represents the constant portion of the source term. With S c fixed, this integral is 


r e r t+At 

S c r"dtdx = S c AtAx p 
J wJ t 1 


(63) 


The fourth integral of Eq. (60) represents the linearly-dependent source term. Here the solution variable i// is 
assumed to be uniform over the cell and fixed at the value at the end of the time step (implicit treatment). The 
resulting integral is 

re rt+At 

S P i/yr"dtdx = S P y/ p AtAx p F (64) 

J w J t 

Combination of Eqs. (60) - (64) results in the following linear equation for the P lh node 

(r//> + j,- T S pAx p)y/ p — R e cip-y/p T R n (tyif/ y + Sp/Ax^, + Up\f/p (65) 

with supporting Eqs. (66) - (70). 
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To support code implementation it is convenient to recast Eq. (65) using notation for the i node as follows 

a,¥i = by i+ 1 + + d t 


(69) 

(70) 

(71) 


The nodal coefficients for the / nodes from N sl ,rf+ 1 to N - 1 are then 



(72) 


&i — cip + R e Ug + R w Ujy *SpAxp 

b, =R e a E 

C i ~ Rn a W 

d t = S c Ax p + a° P i//° P 


Front Boundary Node 

The front boundary condition is applied by imposing a value for the solution variable (temperature or pressure) 
on the front node. This is done by defining the following nodal coefficients for the ;V v „ r / iiodc 


U Nsurf ~ 1 
bNsurf = 0 
^ Nsurf 0 

d Nsurf = V f 


(73) 


where (//, is the value for the solution variable (temperature or pressure) on the front boundary. For pressure, this 
value is simply a specified value. For temperature, this value is iterated until consistency is found between the in- 
depth solution and the boundary condition specification (one of Eqs. (43) and (44)). Details of the iterative surface 
energy balance are given in a subsequent section. 

Back Boundary > Node 

For a specified value of the solution variable (pressure or temperature), the coefficients are simply 


°N ~ 1 

b N =0 

c N = 0 
d N = Vb 


(74) 


where y/b is the value for the solution variable (temperature or pressure) on the back boundary. In addition to a 
specified temperature, general thermal boundary conditions are supported for the back surface, although they are not 
discussed here. The nodal coefficients are found by performing the integration of Eq. (57) over the back element 
(illustrated in Fig. 10). 

Matrix Solution 

The resulting set of N - N sur f + 1 equations of the form of Eq. (71) forms a tridiagonal system as shown in Eq. 
(75). The matrix is solved using the Thomas Algorithm. 



a Nsurf 

b Nsurf 

0 


' C Nsurf +1 

a Nsurf +\ 

_ b Nsurf + 1 

0 

0 

~ C Nsurf +2 

a Nsurf +2 

_ ^ Nsurf +2 

0 

0 



0 


0 

~ C N - 2 

0 

0 


0 

0 

0 

0 



0 

0 

0 

¥ Nsurf 


d Nsurf 


0 

0 

¥ Nsurf +1 


d Nsurf + 1 

0 


0 

¥ Nsurf +2 


d Nsurf +2 


0 

0 


= 


a N- 2 

_ ^N-2 

0 

¥ N- 2 


d N -2 

~ C N- 1 

a N- 1 

-Vi 

¥n- 1 i 


d n- i 

0 

~ C N 

a i V 

¥ N 


d N 


(75) 


Surface Energy Balance 

At each time step, the heat flux conducted into the front surface must be consistent with that determined by the 
front surface boundary condition. The method for reconciling these two terms is described here. 

Heat Flux from the In-Depth Solution 

Solution of the in-depth Energy Equation is accomplished through the solution of the following temperature- 
based form of the system of Eq. (72) 
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Once a surface temperature has been selected, the full matrix system can be solved and the corresponding heat 
flux conducted in-depth q " ond can be determined from the resulting temperature profile. However, solving for the 
heat flux in this manner requires a full matrix solution at each iteration, which is computationally expensive. A 
method is incorporated that avoids the need for a full matrix solution at each iteration. Instead, an expression is 
developed, based on the values of the matrix coefficients, that casts q " cond as a linear function of the surface 
temperature T w associated with the current iteration. This expression then takes the place of the full matrix solution 
in the iteration process. The method is adapted from a similar approach used in the CMA code and is described 
below. 

First, the coefficients for the N sur f node are recast in terms of the unknown surface heat flux q" cond . As a result, the 
definitions of Eq. (73) are temporarily replaced. Integration over the front element yields the following 
modifications for the N sur f coefficients 



( 77 ) 


a Nsurf = a P + R e a E 
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Gaussian elimination is then used to reduce the system to lower-triangular form. The result is 
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For the TV iu ,.^ node, the modified d j coefficient can then be written as 

d Nsurf = RsQcond + $ c AXp + ClpTp + d Nsurf+f Nsurf / ®Nsurf+l 

The equation from the system of Eq. (78) is then 


Rearrangement gives 
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Eq. (83) represents the heat flux conducted into the surface associated with a surface temperature of T w (= T sur f 
Generation of this linear relation requires only one “half-pass” through the matrix to determine the coefficients and 
can be used to represent the in-depth heat flux response to a particular surface temperature during the iteration 


process. 

Surface Heat Flux from Boundary Conditions 

For convenience, a T chem parameter, adapted from previous developers, is defined and incorporated into a table, 
which is referred to as a “Tchem-Table.” The T c ), em parameter is defined, based on the various terms in the surface 
boundary condition, such that the surface heat flux may be written generically as 

C lcond = Pe U e^H (fl r + ^ chem ) + 4 rad,inc ~ £ w <J ^w ( 86 ) 

According to the models of Eqs. (43) - (45), the appropriate definitions are 
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Values for T chem are tabulated versus B' g , B' r . and temperature as illustrated in Table 3. The various terms 
required for T chem are calculated through equilibrium analyses using, for example, the ACE program. The table is 
constructed with predetermined combinations of B' g and B ' c . The temperature, and other thermodynamic conditions, 

are uniquely determined for each set of B' g and B[ . During the solution process, B' g is treated as an independent 
variable with values determined from the in-depth solution at the previous time step, Eq. (15) provides the Darcy 
velocity from which the surface mass flux is determined, and Eq. (46) provides the B' g value. With this in mind, the 

table is organized in sections of constant B' g with varying S c as shown in Table 3. The lines in the table are 
numbered within each B' g section as shown. The remaining terms (temperature, etc.) are uniquely related to the 
particular combination of B' g and Bl c , that is, to the particular line number within a B' g section. The heat flux of Eq. 
(85) can then be thought of as a unique function of the line number of the table within the appropriate B' g section. 



Reconciliation 


Balancing the energy at the front surface is accomplished through the reconciliation of q" ond as calculated by an 
in-depth solution (Eq. (83)) with that calculated according to the specified boundary condition (Eq. (86)) with 
supporting T chem table). That is 

^^Nsurf — P e U e C H {H r +T chem ) (90) 

with T Nsur f set equal to the temperature from the current line of the T chem table. The reconciliation is accomplished 
using a Newton-Raphson scheme to iterate on the line number until the equation is balanced. During this iteration, 
the line number is treated as a real (as opposed to an integer) variable, and linear interpolation against this parameter 
is performed for the dependent variables in the table. Typical convergence occurs at some intermediate (non-integer) 
value for the line number corresponding to an intermediate condition between table entries. Other variables could be 
used as the independent variable instead of the line number, for example temperature or B c . However, there are 
regions within certain types of tables for which one or the other of these parameters may be constant within the 
table. This is the motivation for using the line number as an independent variable. Alternatively, the enthalpy could 
be used. 


Solution Verifications 

The methodologies described above, along with additional capabilities, have been implemented into the ITRAC 
program, and extensive verifications have been performed to ensure that the program provides accurate numerical 
solutions of the governing equations [12]. These verifications have been made by comparing solution results to those 
from analytical solutions, manufactured solutions, and numerical solutions from other accepted programs. In 
addition to the confirmation of solution accuracy, assessment of convergence rates with respect to discretization 
refinement have been investigated. This provides additional verification of proper numerical implementation. An 
exhaustive description of the many verification cases is beyond the scope of this paper. Instead, some selected 
results are summarized. Some comparisons are shown only visually, without detailed descriptions of error norms or 
solution differences. This is done for brevity, with the belief that a visual assessment provides a sufficient level of 
confidence. All verifications shown here that involve pyrolysis have been performed using properties of the 
fictitious TACOT (Theoretical Ablative Composite for Open Testing) material [23]. This material has been defined 
to support evaluation of ablation modeling codes using material properties appropriate for open literature. Six 


verification cases are presented below. These are limited to planer cases highlighting different features of interest. 
Extensive verifications including those in cylindrical and spherical coordinates can be found in the program 
verification reference [12]. 

Case 1 (Pyrolysis Solution) 

Here the focus is on evaluation of pyrolysis (extent-of-reaction) solutions. To evaluate these calculations, a 
constant temperature of 1000 K was imposed on the entire domain of a model, and the extent-of-reaction was 
calculated as a function of time. Pyrolysis kinetics for the charring material were based on the two-component 
model of the TACOT material. Results are compared to those obtained using numerical solution of the governing 
equations (Eqs. (5) and (6)) with Mathematica software [24]. Fig. 11 provides a plot of the two results. The spatial 
domain has no influence on the results here since the temperature versus time history is spatially uniform. 
Accordingly, the results shown correspond to the entire spatial domain. Excellent agreement can be seen, supporting 
verification of pyrolysis solutions. 

Case 2 (Pore Pressure Solution) 

This case provides a comparison of pore pressure solutions against porous media modeling using Fluent® [25], a 
computational fluid dynamics program. Here the pressure was initiated at 100 kPa throughout the entire domain of a 
10 cm slab of a porous material. Pressure conditions on the front and back surfaces were prescribed at 200 and 100 
kPa, and the transient response was modeled. Material porosity and permeability were fixed at constant values of 0. 1 
and 1 x 10" L ’ m 2 , respectively. Viscosity and molecular weight for the flowing fluid were based on fixed values for 
air at 300 K. Results of the calculated pore pressure response are shown in Fig. 12 along with Fluent results. The 
plots show pressure profiles at various times. Excellent agreement is seen, supporting verification of the pore 
pressure solutions. 

Case 3 (Temperature with Moving Boundary Solution) 

Here comparisons are made to an analytical solution. A constant surface temperature T s of 2000 K and a constant 
erosion rate s of 1.27 mm/s were applied to the front surface. A thickness of 0.6 m was modeled, sufficient to ensure 
semi-infinite behavior. The entire domain was initialized at a temperature 7] of 300 K, and constant thermal 
properties were used. The simulation was run for a model time of 400 s, and discretization was based on 2000 
elements of equal thickness and a time step size of 1 s. As time progresses in this type of model, the temperatures 


ahead of the moving surface approach a quasi-steady profde for which an analytical solution can be derived. 
Carslaw and Jaeger [26] provide the analytical solution, which can be written as 

T = T i + T s e~ ix ' la (91) 

where x r is the depth relative to the front surface and a is the thermal diffusivity. Fig. 13 provides plots of the 
solutions, in the form of temperature profdes, versus those based on Eq. (91). The times compared are sufficient for 
quasi-steady conditions to be reached. The comparison shows excellent agreement for this thermal solution 
including the effect of a changing domain. 

Case 4 (Pyrolysis, Pore Pressure, and Temperature Solutions) 

This case considers coupled solutions for in-depth heat transfer, material pyrolysis, and pore pressure. The 
problem definition is that used as an initial test case for code comparisons as part of the recent Ablation Modeling 
Workshop [27]. A surface temperature of 1664 K was applied to one side of a 5 cm domain, while the other side 
remained adiabatic. The entire domain was initialized at a temperature of 298 K, and the simulation was run for a 60 
second duration. Properties for the TACOT material were used. No analytical solution is available for this type of 
test case. Instead, results are compared with those from the familiar CMA code to ensure that solutions are 
consistent with a commonly accepted approach. The solutions were obtained using 1000 elements of equal thickness 
and time steps of 0.01 s. CMA [16] solutions were run with 600 elements and a maximum time step of 0.1 s. Both 
models were evaluated for discretization convergence. Results are shown in Fig. 14 and Fig. 15. The first figure 
shows temperature comparisons at six depths relative to the front surface. The second figure shows depths of the 
pyrolysis and char fronts, defined as 2% and 98% charred, respectively. Excellent agreement is seen between the 
two codes. This verifies consistency with an accepted solution method for charring materials including the effects of 
heat transfer, material pyrolysis, and pyrolysis gas flow. 

Case 5 (Pyrolysis, Pore Pressure, Temperature, and Thermochemical Erosion) 

This is another case from the Ablation Modeling Workshop [28]. In addition to material charring, this case 
includes the effects of surface thermochemistry and erosion. A 5 cm thick sample of the TACOT material was 
heated on one side by thermochemical convection from air flow, while the other side remained adiabatic. The heat 
transfer coefficient was ramped linearly from an initial value of zero to its full value at 0.1 s, then held constant 
through 60 s. The value of the enthalpy-based heat transfer coefficient was 0.3 kg/m 2 -s, and the corresponding 



recovery enthalpy was 2.5 x 10 7 J/kg. The entire domain was initialized at a temperature and pressure of 300 K and 
101.325 kPa. Supporting thermochemistry (B-prime) tables were created using a modified ACE program [17] with 
air chemistry for the adjacent boundary flow. Mole fractions of 0.21 and 0.79 were used for oxygen and nitrogen, 
respectively, and the unity Le model was used. The insulator and pyrolysis gas chemistry were defined based on that 
for the TACOT material. Verifications were made against solutions from the CMA program. Discretizations were 
based on 1000 elements and 0.01 s time steps. The CMA model used 600 elements and a 0.1 s maximum time step. 
Both models were evaluated to ensure discretization convergence. Results showing temperatures at various depths 
from the original surface location are shown in Fig. 16. Results for the eroded surface, the pyrolysis front (2% 
charred), and the char front (98% charred) are shown in Fig. 17. Excellent agreement is seen between the two 
models verifying consistent solutions with accepted methods for a charring material with chemical surface ablation. 


Case 6 (Discretization Convergence Rates) 

The rate of solution convergence with respect to discretization refinement is assessed in this verification case, 
which considers heat transfer with constant properties and a moving surface. The convergence rate is evaluated 
based on the behavior of the L 2 norm of the error with respect to an analytical solution. Specifically, derivatives of 
the Li norm with respect to spatial and temporal discretization sizes are evaluated to confirm that the error norm 
converges toward zero at expected rates. The method of manufactured solutions [29], [30] is used in the assessment. 
Method of Manufactured Solutions 

An analytical solution is obtained using the method of manufactured solutions. With this approach, a solution 
rather than a problem statement is chosen. The corresponding problem statement is then derived using a variable 
source term in the governing equation. In this manner, a problem statement with its corresponding analytical 
solution is obtained. The method is applied here including the effect of a moving boundary. 

For planar conditions with constant properties Eq. (59) can be rewritten as 


c^L.k^s, 

dt dx 2 c 



(93) 



Eq. (93) represents the equation solved in this verification case with q calculated as a function of space and time 
according to the manufactured solution described below. 

The changing spatial domain is illustrated in Fig. 18, which shows the original location of the front (ablating) 
surface at x = 0 . the current location of the front surface at x = x v , and the fixed location of the back surface at 
x = x b . The overall length of the spatial domain L, which changes with time due to the movement of x. , is 

L = x b -x s (94) 

Nondimensional spatial coordinates t// and <f> are defined for convenience. These are illustrated in Fig. 19 and are 
defined as 
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and 
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In addition, the time domain is nondimensionalized with respect to the total (final) simulation time t f as follows 


9 = — 


(97) 


A manufactured temperature solution T m is then defined using these nondimensional coordinates as 

T m =a ¥ b f9 d (98) 

For an instantaneous ablation rate of s , the temporal and second spatial derivatives, with respect to the x and t 


coordinates, are 
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With constant properties, substitution into Eq. (93) provides the required source term as 



(101) 
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where a is the thermal diffusivity. 

For this case, the exponents b, c, and d in Eq. (98) were chosen to be 2, 2, and 1 respectively, and the coefficient 
a was 1 X 10 s . This provides an exact solution with a constant temperature of zero at the endpoints of the domain 
(even as ablation progresses) and a quartic trace on the interior that increases linearly with time. The original length 
L of the domain is 0.1 m, the ablation rate s is held constant at 5 X 10 -4 m/sec, and the thermal diffusivity a is 
1 x 10~ 5 m 2 /sec. 

Error Convergence 

Error in the numerical solution at any point in time and space is defined as the difference between the 
numerical solution T and the exact solution T m from Eq. (98), that is 

e(t,x) = T(t,x)-T m (t,x ) (1Q2) 

At any point in time, the overall spatial error norm is calculated as 

(103) 

This is integrated over the temporal domain to provide the total error as follows 

(104) 

Solution convergence is assessed by evaluating the total error versus discretization levels. Derivation of the 
discretization equations employed a zeroth-order (constant) temporal representation and a first-order (linear) spatial 
representation. As a result, errors associated with discretization should be expected to follow 

llell =aAt + bAx 2 (105) 

II iil 2 

In regions where temporal error dominates (highly refined Ar) this leads to 

logle^ = log a + log At (106) 

When spatial error dominates (highly refined At) the expected behavior is 

logllell = logZ> + 21ogAx (107) 

II IIL 2 

Differentiation provides the following when temporal errors dominate 





51 °g|HL 2 1 

d log At 

{first order) 

(108) 

When spatial errors dominate the relation is 



ai °g|HL 2 , 

d log Ar 

(s “" ) 

(109) 


The expected rates of convergence for temporal and spatial discretization are then 1 and 2, respectively. Error norms 
have been calculated with various discretization levels to assess the actual convergence rates. Gaussian quadrature 
was used for the integration of the error norm. 

Results 

The numerical verification of the variable grid method was accomplished using an array of temporal and spatial 
discretization values. These values were selected to span a range for the specified source term, derived as part of the 
method of manufactured solutions. Fig. 20 and Fig. 21 illustrate results for coarse and fine spatial discretization, 
respectively. For the coarse model of Fig. 20, initial spatial discretization included 10 elements with uniform nodal 
$pa cing of 0.01 m while temporal discretization used equal times steps of 1.0 sec. For the refined model of Fig. 21, 
spatial discretization included 80 elements with 0.00125 m nodal spacing and 0.1 sec time steps. As expected, 
improvement is seen with an increased level of discretization. To assess the convergence rates, error norms were 
calculated with various levels of refinement. Table 4 shows results where temporal error dominates. Time steps here 
range, in factors of 10, from 10 to 0.001 sec. Logs of the time step and error are listed in the table along with the 
numerical derivative representing the convergence rate. Similar results for spatial refinement are given in Table 5; 
here, spatial effects were dominant. The numerical derivatives listed in the tables show the expected convergence 
rates of 1 and 2. The accuracy illustrated in Fig. 21, along with the proper convergence behavior, provide evidence 
for accuracy and proper code implementation, including the effects of a moving boundary. 

Solution Example 

A unique advantage of the present method is the ease with which surface movement can be controlled. This 
supports, for example, the modeling of materials with char layers that fail under certain structural loading 
conditions. To help investigate related effects, several mechanical erosion submodels are included in the ITRAC 
program. These include a limiting thickness for the char layer, an instantaneous spallation of the char layer after 
reaching a spallation criterion, and a mechanical erosion rate that augments the chemical rate. In any of these 



models, the ablation augmentation is accounted for by simply moving the surface to the appropriate location using 
the grid modification scheme previously described. 

A common need for mechanical erosion modeling is related to elastomeric insulators such as those used 
internally in solid rocket motors. Char characteristics for these materials are notoriously weak and difficult to 
characterize. As a result, successful models are often calibrated to available test data, and the capabilities of the 
present solution scheme allow for these investigations to be performed in a meaningful way. As a demonstration, 
consider the behavior of an elastomeric insulator in a subscale solid rocket test motor where the material is known to 
exhibit a pyrolysis depth of 7.6 mm at a particular station with an exposure time of 48 sec. The material behavior 
and boundary conditions are hypothesized, but are typical of actual material behavior. Predictions using 
thermochemical surface ablation give the transient response shown in Fig. 22. The figure shows the modeled depths 
for the surface, the char line, and the pyrolysis line, where the char and pyrolysis lines correspond to 98% and 2% 
extent-of-reaction, respectively. Here the surface has moved, due to thermochemical surface ablation, to a depth 
around 2.7 mm and the overall pyrolysis depth is just under 4.0 mm. It is speculated that the deeper pyrolysis depth 
of 7.6 mm is a result of mechanical material failure and a model is constructed that forces a spallation event to 
periodically occur. In this case the spallation is based on a limiting char layer thickness; whenever the char layer 
reaches a specified thickness, a spallation event removes any material with an extent-of-reaction higher than 0.5. 
This is accomplished within the model by simply moving the surface to the appropriate depth when the criterion is 
met. Results including the spallation model are shown in Fig. 23. The spallation event can be seen in the program 
output, and the pyrolysis depth is augmented to a value near the 7.6 mm observation. Surface temperatures 
corresponding to the two models are shown in Fig. 24. This simple example illustrates the usefulness of the 
modeling approach for modeling and investigating complex ablation phenomena including effects of mechanical 
erosion. 


Discussion 

•The numerical method described herein has been used to create the Insulation Thermal Response and Ablation 
Code, or ITRAC. Only a fraction of the capabilities of the ITRAC program have been presented in this paper, and 
the verifications presented have been limited to those highlighting general solution accuracy and comparisons with 
the long-accepted CMA program. Those comparisons have focused on problems for which agreement with CMA 


solutions should be expected. Details of additional solution capabilities and extensive verifications can be found in 
complete theory, user, and verification manuals [10] - [12]. The ITRAC program provides a modern alternative 
lacking the numerical issues commonly encountered with the CMA program; these include unexpected solution 
instabilities and failure to converge with discretization refinement under certain conditions. Other advantages of the 
ITRAC program include various modeling options, some of which have been described here. For brevity, many of 
the unique features have been omitted from this discussion entirely, but they include various options for mechanical 
erosion, hill pore pressure solutions for applications where related mechanical loading is of importance, multi- 
component pyrolysis analysis with no limits on the number of components, alternative surface ablation models 
including heat-of-ablation techniques, temperature-based boundary conditions, multiple options for defining 
property dependence, an option for implicit solution coupling, internal diffusion of absorbed constituents such as 
moisture, and capabilities for user-defined thermal and erosion boundary conditions based on other solution 
parameters. The program is extensively used within the ablation modeling community and consistently provides 
stable solutions even with conditions of extreme heat flux and ablation conditions. For complete descriptions of the 
program, the reader is referred to the ITRAC manuals. 

Summary and Conclusion 

Ablation heat transfer modeling equations have been presented. These mathematical models include the effects 
of heat transfer, material pyrolysis, internal permeation and heat exchange, and thermochemical surface ablation. 
Detailed derivations have been provided so that modeling assumptions are clear, and the models have been put into a 
form that supports coupled numerical solution. A one-dimensional numerical solution scheme has also been 
presented. The scheme is based on a control-volume formulation with a variable grid to account for the effects of 
surface movement. The solution method has been implemented into a new computer program and accuracy of the 
method has been verified. The variable grid method used in the program allows for easy implementation of surface 
movement models associated with phenomena such as mechanical erosion. A simple case was presented 
demonstrating this capability through the application of a spallation model used to simulate periodic mechanical 
char removal of an ablative. This general modeling approach is expected to provide a foundation for continued 
development and improvements in the important area of ablation heat transfer. 
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Fig. 5 Control-surface for surface energy balance. 


Ablating surface, S 
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Fig. 7 Control-volume modifications at the ablating surface. 
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Fig. 14 Temperature results for verification case 4. 
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Fig. 16 Temperature results for verification case 5. 
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Fig. 19 Nondimensional coordinates. 
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Fig. 20 Results for verification case 6 with coarse discretization. 
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Table 1 Generic parameter definitions for the energy and mass/momentum equations. 





Table 2 Curvature correction factor definitions. 


Curvature 

F r 

Planar geometry ( li = 0 ) 

1 

Cylindrical geometry ( n — 1 ) 

r e + r w 
2 

Spherical geometry ( n = 2 ) 

2 2 
r e +r e r w +r w 

2 




Table 3 T chcm table format. 






Table 4 Error norms with temporal dominance (Ax = 0.00002 m). 


At 

log At 

log E 

Alog(E)/Alog(At) 

10 

1 

2.925341 


1 

0 

1.982851 

0.942490 

0.1 

-1 

0.988789 

0.994062 

0.01 

-2 

- 0.010576 

0.999365 

0.001 

-3 

- 1.010097 

0.999521 





Table 5 Error norms with spatial dominance (At = 0.001 s). 


Ay 

log Ay 

log E 

Alog(E)/Alog(Av) 

0.01 

-2.000000 

2.498373 


0.005 

-2.301030 

1.872606 

2.078753 

0.0025 

-2.602060 

1.251767 

2.062383 

0.00125 

-2.903090 

0.638797 

2.036242 

0.000625 

-3.204120 

0.033858 

2.009564 




